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We study proximity effects in clean nanoscale superconductor-normal metal-superconductor (S|N|S) 
graphene heterostructures using a self-consistent numerical solution to the continuum Dirac Bogoliubov-de 
Gennes (DBdG) equations. We obtain results for the pair amplitude and the local density of states (DOS), as 
a function of doping and of the geometrical parameters determining the width of the structures. The supercon- 
ducting correlations are found to penetrate the normal graphene layers even when there is extreme mismatch in 
the normal and superconducting doping levels, where specular Andreev reflection dominates. The local DOS 
exhibits peculiar features, which we discuss, arising from the Dirac cone dispersion relation and from the in- 
terplay between the superconducting and Thouless energy scales. The corresponding characteristic energies 
emerge in the form of resonant peaks in the local DOS, that depend strongly on the doping level, as does the 
energy gap, which declines sharply as the relative difference in doping between the S and N regions is reduced. 
We also linearize the DBdG equations and develop an essentially analytical method that determines the critical 
temperature Tc of an S|N|S nanostructure self-consistently. We find that for S regions that occupy a fraction of 
the coherence length, Tc can undergo substantial variations as a function of the relative doping. At finite tem- 
peratures and by manipulating the doping levels, the self consistent pair amplitudes reveal dramatic transitions 
between a superconducting and resistive normal state of the structure. Such behavior suggests the possibility 
of using the proposed system as a carbon-based superconducting switch, turning superconductivity on or off by 
tuning the relative doping levels. 
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1. INTRODUCTION. 



The successful development of methods to create large 
samples of graphene,'-^ has been followed by recent efforts 
to exploit its high electron mobility-' "* and the peculiar band 
structure^ associated with its two dimensionality. A number 
of graphene based devices have been subsequently proposed*", 
including field effect transistors, quantum information storage 
systems, optoelectronic devices, and nanoscale superconduct- 
ing systems. In particular, the observation of superconduc- 
tivity in graphene,^"' ' either through doping or by means of 
superconducting contacts, has fueled research activity involv- 
ing proximity effects in normal (N) and superconductor (S) 
graphene regions that are in close electrical contact. Indeed, 
the presence of superconducting correlations in graphene is 
remarkable considering that undoped graphene in isolation 
is inherently nonsuperconducting even at low temperatures. 
Striking evidence of the peculiarities of superconductivity in 
graphene was the observation of a Josephson supercurrent in- 
duced by two superconducting electrodes in close contact with 
the graphene.** The tunneling conductance of a junction con- 
sisting of an insulating barrier between graphene and a super- 
conductor should exhibit oscillations'-' as a function of bar- 
rier strength, surprisingly peaking at finite values. These un- 
expected effects arise in large part from the hexagonal sym- 
metry of graphene, which generates a relativistic-like band 
structure^ near six points on the Fermi surface-the so called 
Dirac points. The low energy dispersion near these points is 
linear, and subsequently, the quasiparticles are governed by a 
two dimensional massless Dirac-like equation. 



Superconducting proximity effects in conventional het- 
erostructures consisting of a normal metal and superconductor 
have been known for a very long time. If the superconductor 
is coupled to a graphene sheet, where "Dirac quasiparticles" 
are confined to a 2-D plane, the leakage of superconductiv- 
ity into graphene should exhibit novel behavior Thus, study- 
ing superconducting proximity effects in graphene requires 
a careful and accurate determination of the pair correlations 
throughout the entire system. These are characterized by the 
pair potential A(r), and the pair amplitude, F{r). A proper 
delineation of the associated proximity effects can only be 
achieved through a self consistent calculation of A(r), ensur- 
ing that the system's lowest free energy state is found. The 
resultant self consistent state generally possesses nontrivial 
spatial inhomogeneity that can have important consequences 
for quasiparticle bound states, interface bound states'^ and 
potential supercurrent flow. It is not surprising then that the 
frequently used step-function model for A(r), while satisfac- 
tory for length scales much longer than the superconducting 
coherence length, ^q, can lead to erroneous results for small 
graphene structures where quantum scale oscillations play a 
role. For example, superconductor widths that are on the 
same order as give rise to self consistent pair potentials that 
can vary over a significant fraction of the total sample width. 
Moreover, self consistency is crucial at finite temperatures, 
where the superconducting correlations can have substantial 
decay near the interfaces. 

The usual superconducting proximity effect is governed by 
the mechanism of Andreev reflection. This is the process 
where at the interface, an electron with energy below the su- 
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perconducting energy gap is retro-reflected as a hole, transmit- 
ting a Cooper pair into the superconductor In graphene, the 
effectiveness of the Andreev process depends in part on the 
relative doping in the S and N regions. For electron doping 
the Fermi level is shifted upwards, while for hole doping, it is 
shifted downwards relative to the Dirac point. If the normal 
graphene layer is weakly doped, so that its Fermi level, n^, 
is in absolute value much smaller than that in the S region 
(|mjv/msI ^ 1), ipecM/ar Andreev reflection becomes impor- 
tant. In this process, the electron and hole belong to differ- 
ent bands. Thus, despite large Fermi wavevector mismatch, 
superconducting correlations can penetrate into the normal 
graphene region. If on the other hand, fiN/fJ-s ^ both 
the conventional and specular Andreev reflection processes 
are suppressed and normal scattering drives the quasiparticle 
trajectories. The doping level clearly then has important con- 
sequences for any thermodynamic and transport properties in- 
volving superconducting graphene nanojunctions. 

Besides doping effects, there are geometrical issues to con- 
tend with in finite S|N nanojunctions: the electronic structure 
of confined graphene can lead to a strong size dependence.'^ 
Depending on the widths of the normal graphene and super- 
conducting regions, there are various energy scales that can be 
difficult to disentangle. If a S|N|S heterostructure has a thin 
middle channel of width much smaller than (the super- 
conducting coherence length), the relevant low energy scale is 
the usual energy gap, Aq. This holds then for superconducting 
graphene'^ which behaves, in this case and in this respect, the 
same way as conventional three-dimensional materials. '^'^"^ 
Short structures also result in critical currents that can devi- 
ate from the simple harmonic form.^' For wide middle layers, 
with cIn ^ £,Q, the Thouless energy, Et = Vf/d^, (vf is 
the Fermi velocity) emerges as an important energy scale. 
The Thouless energy in clean systems gives rise to geometry 
dependent quantum phenomena that arise from the cumula- 
tive phase coherent effect of propagation and reflections from 
the structure boundaries. This energy scale interacts then with 
the Ao scale: for large normal graphene widths, Et can be 
smaller than Aq, and the energy spectrum possess a Thou- 
less gap for quasiparticles with energies less than a charac- 
teristic energy of order Et- When Et is of the same or- 
der as Aq, identifying the origin of spectral anomalies can 
be difficult. For long S|N|S heterostructures, the resultant 
peaks in quasiparticle spectra would be (assuming a non self- 
consistent, step function form for A(r)) located at energies 
proportional to integer multiples of Et-^^ Self consistency can 
modify this result however as the pair potential deviates sub- 
stantially from a simple step function model. Moreover, when 
the doping amount changes, this picture becomes complicated 
by changes in quasiparticle bound states and density of states 
(DOS) due to the shifting of the Fermi level. 

In this paper we use a fully self consistent framework to 
calculate the energy spectrum and pair amplitude in S|N|S 
graphene heterostructures. In addition to the self consis- 
tent pair amplitude, our accurate numerical diagonalization 
method allows us to investigate two important quantities that 
can be measured experimentally and further clarify proxim- 
ity effects in graphene. The first is the local density of states 



(DOS), which can be measured directly with a scanning tun- 
neling microscope. The vanishing of the DOS at the Fermi 
energy in undoped graphene results in varying subgap bound 
states and minigaps^-' associated with the interplay between 
Et and Aq. To effectively characterize the local electronic 
properties, we examine the DOS, in both the S and N regions. 
The energy spectra reveal conditions for fully gapped and gap- 
less states in graphene S|N|S junctions. For a given doping 
level, the gap width and magnitude are shown to diminish as 
dN increases. The greatest variations are found to occur when 
Et < Aq. 

The second experimentally observable quantity of interest 
is the critical temperature, Tc- We study how Tc varies as a 
function of doping levels and of the geometrical parameters. 
Our self consistent calculations find a nontrivial variation in 
Tc as a function of the relative doping levels, hn/ l^s- The 
sensitivity of Tc on the Fermi shifts depends strongly on the 
width of the outer S regions: very thin superconductors with 
ds < £.0 reveal the most drastic changes in Tc for small in- 
crements in doping. We show that, for particular ranges of 
ds, c?Af, and temperature, a S|N|S nanostructure can act as a 
type of switch that transitions between a superconducting and 
resistive normal state as the ratio /iAr//^5 is varied, something 
that might be done by using a modulated in-plane external 
electric field ''^"*. 

Despite the importance of self-consistency, the only pre- 
vious self-consistent works addressing proximity effects in 
S|N|S structures are Josephson junction studies based on an 
extended Hubbard model with the superconductivity arising 
either from doping" or from external contacts. ' It is the 
aim of this paper to present a method that complements tight 
binding approaches and provides a suitable description for 
the critical temperature and applicable for a wide range of 
geometrical and coherence lengths. We achieve this goal 
by numerically solving the microscopic Dirac Bogoliubov-de 
Gennes (DBdG) equations self consistently in the continuum 
regime. By retaining atomic length scales in the calculations, 
we can accurately represent the important geometrical effects 
inherent to finite sized junctions. The DBdG equations are 
ideal for inhomogeneous S|N|S heterostructures since they 
give directly the quasiparticle amplitudes and energies that 
characterize proximity effects. They are also appropriate for 
clean systems such as graphene which has high electron mo- 
bility. To investigate the potential of our S|N|S system as a 
superconducting graphene switch, we determine the critical 
temperature self consistently. This is accomplished by taking 
the full DBdG equations, and linearizing them via standard 
perturbative techniques. We then arrive at an essentially an- 
alytical method that determines the critical temperature as a 
function of geometrical parameters and doping levels. 



2. METHOD 

The geometry we study consists of a graphene sheet in- 
finite in one direction (that of the y axis) and comprised of 
two doped strips of superconducting material, each of width 
ds, separated by a normal region of width rf^v- We consider 
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the pairing in the S regions to be conventional s-wave. The 
methods we use to self-consistently diagonaHze the mean field 
single-band Hamiltonian are extensions of those previously 
employed^^^' to study proximity effects in ordinary three di- 
mensional materials, but important changes have to be made 
to take into account the reduced dimensionality and the pecu- 
liar band structure of graphene. These changes are the focus 
of the discussion below. 

Our starting point in this case is the Dirac Bogoliubov- 
de Gennes (DBdG) equations which govern the quasiparticle 
spectrum of graphene. In the absence of magnetic effects, the 
DBdG equations for the two valleys K{+) and A''(— ) are'*: 
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The Dirac Hamiltonian, 'H^'^, is given compactly hy Ti^ = 
Vf{(JxPx ± '^yPy)^ in which (Ji are the 2x2 Pauli matri- 
ces acting in sublattice space, / is the identity matrix, w/ is 
the (energy independent) Fermi velocity in graphene, and /i 
is the chemical potential. This quantity vanishes in the un- 
doped case but not in the presence of doping. Here we will 
take ^{x) to be a piecewise constant: a fixed positive number 
Us in the S region and a variable value ji^ in the regions. 
We consider the case of relatively large doping in the S re- 
gions, so that we can assume smooth interfaces.'*'-* We will 
take these interfaces in the direction of constant x. We have 



defined *+„ = 



'■B,K'l 



.... K,K,«S,A')'^ and ^ ivli,,,v"sj,,y- The 
A, B labels denote the two sublattices that arise from the hon- 
eycomb lattice structure. 

The notation in the Hamiltonian implies that the Pauli ma- 
trices act on the pseudospin of the quasiparticles, mapping the 
usual spin into the projection of the wavefunction onto sub- 
lattice A or B. Since the valleys are degenerate {K and —K' 
are equivalent), we need only solve for either or . As- 
suming the first choice, we define the four component vector 

^ri = (^u,n'^u,n)- 

The pair potential A couples electrons in a given valley with 
the hole excitations in the other valley. In terms of the wave- 
functions and energies obtained from Eq. (1), this coupling 
leads to the self-consistency condition, 

A(^) = f E K.kvXk' + n%Kvi:K] taiih(^) , (2) 

n 

where the superconducting coupling parameter, g, is a positive 
constant in the intrinsically superconducting regions and zero 
elsewhere. The sum is over all energy eigenstates within the 
Brillouin zone whose energy, referred to jis, is smaller than 
or equal to a characteristic energy cutoff, It is to be inter- 
preted as — > 1/ (27r) J dky where ky is the transverse 
momentum, and q a longitudinal index. The singlet pairing 
only occurs from opposite valleys, to maintain time-reversal 
symmetry'^''*. 

The experimentally important local DOS, N{x, e), is given 



by 

N{x,e) 



n,a,/3 



(3) 



where a equals Aox B, (3 can be either K or K' , and /' is the 
derivative of the Fermi function. One can integrate N{x, e) 
over any suitable range of x to obtain the average DOS in a 
certain region. 

We now take advantage of the translational invariance along 
y by writing 5'„(a;,?/) = e*'^"^$„(x). Introducing the nota- 
tion <l>^(x) = (s„(2;),t„(a;), w„(a;), z„(a;))-^ where the func- 
tions in the parenthesis correspond to valley and sublattice in- 
dices in the same order as for the previously defined 5'„, we 
can rewrite the DBdG equation Eq. (1) as. 
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where we define ti± = ivf{dx ^ ky), (we use h = 1 and 
kg ~ 1 throughout this paper). 

Next we expand the quasiparticle wavefunctions via, 
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'^nix) ^'^Cn,q(pq{x), 



(5) 
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are the expansion 



where the c„ 

coefficients, and the (j)q (x) is a set of N basis functions, where 
must be sufficiently large. ^^-^^ We take, consistent with 
the boundary conditions, (t>q{x) = ^y2/dsm{kqx) in which 
d = 2ds + dN, and kq = qir/dis the quantized wavenumber 
These basis functions are not eigenstates of the normal Hamil- 
tonian. Therefore there are tt^ off diagonal terms. These in- 
troduce some computational challenges that result in a larger 
value of N being required than in the three dimensional cases. 
Considering then each row of Eq. (4), we perform the fol- 
lowing somewhat lengthy but elementary steps. First, we in- 
sert the expansion Eq. (5) into Eq. (4). Next, we multiply 
each term by (pqi and integrate the variable x over the range 
< X < d, taking into account properly the stepwise x de- 
pendence of fi. Finally, we choose i^is as our unit of energy 
(recall that we are dealing with a strongly doped system so 
this quantity does not vanish) and divide through by 

Taking the same steps for the rest of the matrix, we end up 
with the following 4N x AN matrix equation: 
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and the vector, a„, contains the expansion coefficients, 

Wn,l,--- ,Wn,N,Zn,l, - ■ ■ , Zu^n)'^ ■ (8) 

Here, I and O are unit and zero matrices of rank N re- 
spectively. Consistent with our choice of energy units, we 
now define tilded dimensionless energies Jl^ = Hn I f^s 
?„ = €n/ ^s- We choose also to measure our wavectors in 
units of kps defined by the relation kps = ^s/vj, thus e.g. 
ky = ky/kps and then have for the remaining elements in 
Eq. (7): 

"4<j,<j' = ICq+q'{ds) - ICq^q'{ds) - JlN^Cq-q' {ds + d^) 
- ICq+qi{ds + d^) - ICq^q'{ds) + JCq+q>{ds)] 
+ ICq-q'{ds + dN) - ICq+q'{ds + dN), Q ^ q' , (9) 
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ds ~ r 



■^2q{ds) 



-K,2q{ds + dN)]- - lC2q{ds + d^), (10) 



2iq'q ( -1 + 
kpsd y - q'^ 



I, q^q\ 
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and Bq^q = 0. Here, /C„(a;) = sin(n7ra;/(i)/(n7r). Finally; 



new A(.t) from Eq.(13) and iterate this process until the rela- 
tive difference between successive A(a;) is less than 10^^. 

We also determine the critical temperature, Tc, semi- 
analytically as a function of the geometrical and doping pa- 
rameters. To find Tc, the self-consistency equation can be 
linearized^' near the transition, leading to the form 
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"^2(7 Aq , 
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where the A^ are expansion coefficients of the pair potential 
(Eqs. (2) and (13)) in our basis and the Jiq are the appropriate 
matrix elements with respect to the same basis, as obtained 
from the linearization procedure. These matrix elements can 
be written as Jiq = (J*^ + Jiq)/2, where. 
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(17) 
(18) 



Here 7 = \/ {2'iT'^kFsd), with A the dimensionless supercon- 
ducting coupling constant introduced above. The eigenen- 
ergies, en^^^'^, are the unperturbed particle (hole) energies 
(found by setting A = in Eq. (4)), and Nd denotes that 
the sum is cut off at energies beyond the ojc frequency. We 
also have. 
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where Aq is the order parameter in bulk S material at T = 0. 
The self-consistency relation Eq. (2) can now be written as. 



A(a;)/Ao = 4A 
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where the correlation factor, ICqpr, is written. 
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+ 'tn,qZ*,^q,) sm{kqx) sm{k'^x) tanh(^|^^ , 

where^*^ ^0 — Vf/Ao, kc is the ky cutoff corresponding to 
those states specified below Eq. (2), and A is the dimensionless 
coupling constant which we define as A = giJ,s/2TTVj. Since 

the DOS for bulk S material in its normal state is Noi^-) = unity^''-^'* corresponds to Tc. This linearization route is much 
2e/{Trvj), we have A = .gA^o(/"s)/47r. 

To perform our calculations, we must solve Eq. (6) to- 
gether with the self consistency condition Eq. (13). When 
ojc = Wc/ /is satisfies tj^ < 1, we find for the bulk case. 



Here we define &{z) to be unity in the superconducting re- 
gions and vanish in the normal ones. The determination of 
Tc involves calculating the eigenvalues of matrix Jiq for each 
temperature value in the range of interest, and the highest tem- 
perature for which the largest eigenvalue of the matrix Jiq is 
,33,34 

more efficient at calculating Tc than solving the full DBdG 
equations near Tc, which can often be very difficult due to 
the large number of iterations involved in the self consistency 
process. 



A = arcsinh(a;c/Ao), 
while if cDc > 1, then-*^ 



(14) 



A-^ = yA2+w2_ ^A2 + l + arcsinh(l/Ao). (15) 

where Ao = Aq/hs- To achieve self consistence, we start 
with an initial guess for A (a;) and once all of the eigenfunc- 
tions and eigenenergies have been determined, we calculate a 



3. RESULTS AND DISCUSSION 

In this section, we present our self consistent results for 
the pair amplitude, local DOS, and critical temperature, for a 
broad range of widths and relative doping levels. We con- 
sider the S regions to be electron-doped, corresponding to 
/i5 > 0, while the normal graphene can be either electron- 
doped {pLM > 0) or hole-doped (/ijv < 0). All lengths are 
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FIG. 1: (Color online). Normalized pair amplitude vs X (see text) at T = for an S|N|S heterostructure. The S portions have a dimensionless 
width (see text) of Ds = 150 each. For clarity, the outermost parts of the sample are not shown. Each panel corresponds to a different 
normalized normal graphene width of (a.) Dn ~ 20, (b) Djv = 50, (c) Dn ~ 100, and (d) Djv — 300. For each case, four curves 
representing different doping levels are shown. From top to bottom in the central N region (green, blue, black, red), these curves correspond 
to JIn = 0.5, 0.2, 0, and 10. 



scaled in units of the Fermi wave vector kps, and we define 
the relative dimensionless coordinate X = kpsix — d/2), so 
that X = is at the center of the structure. When consider- 
ing thermal effects, all quantities involving the temperature, 
T, are scaled by Tq, the transition temperature for the bulk su- 
perconducting material. Our input parameters are, besides the 
geometrical lengths, which are given in dimensionless form as 
Ds = kpsds and Dn = kpsdN, the value of ptjv, and that 
of the dimensionless coherence length So = ^fs^Co = 100. 
From the latter, and using a fixed value of lJc = 0.04 (for 
U!c < 1, results are only weakly sensitive to uJc) we obtain A 
viaEq. (14). 

We will first consider our results for the self-consistent 
normalized Cooper pair amplitude F{x) = (l/A)A(a:;)/Ao, 
which reveals the superconducting correlations throughout the 
entire S|N|S system. Some of our results for F{X) as a func- 
tion of dimensionless distance X are shown in Fig. 1 and 
Fig. 2. In Fig. 1, F{X) is shown, in each panel, for a different 
value of Dn, and all at the same Ds = 150. Several (all pos- 
itive) values of the relative doping parameter Jl^, correspond- 
ing to electron doping, are shown in this figure: the undoped 
N case is also shown for comparison. We see that the proxim- 



ity effect depends strongly on the relative doping JIn via the 
mismatch it reflects (when this quantity is unity, there is no 
mismatch). Moderate doping in the N region allows for very 
appreciable pairing coiTelations in the normal graphene, even 
when Dpf > Sq (right bottom panel), as one can see e.g. in the 
JiN — 0.5 ((green) highest curves at X — 0) results shown, 
while depletion of F{X) in the S regions extends in this case 
to distances longer than the correlation length. On the other 
hand, when the mismatch in the Fermi shifts is extreme (as in 
the /iAT = ((black) solid curve, third from top at X = 0) and 
/iw = 10 ((red) dashed curves) cases shown, we see that the 
proximity effect is much weaker The (blue) curves, second 
highest at X ^ corresponding to /ijv = 0.2, show interme- 
diate behavior Further examination of the results in this fig- 
ure for large mismatch CjlN — 0) reveal that specular Andreev 
reflection allows somewhat more readily for the penetration 
of coiTelations in the normal graphene than large mismatch in 
the opposite direction: comparing the /Iat = to the ptjv = 10 
results, the self consistent state in the later case shows less su- 
perconducting correlations in the normal graphene than in the 
former case, and a correspondingly smaller depletion in the 
S layers. The sequence of panels, moreover, illustrates that 
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FIG. 2: (Color online). Normalized pair amplitude versus X, as in Fig. 1. The dimensionless geometrical quantities Dn and Ds are now 
Dm = 300 and Ds ~ 300, so that the S|N boundaries are at \X\ = 150. The top set of panels depicts most of the S|N|S structure and the 
bottom set of panels are magnifications of part of the N region near the interface for the panel above. Panels (a) and (c) are for electron doped 
N with curves shown, from top to bottom al X = 0, for /Ijv — 0.5, 0.2, and 0.1 (red, green and blue), while panels (b) and (d) are for hole 
doped N with ptjv — —0.1, —0.2, and —0.5 with a similar scheme. In both cases results for ptjv — (black) lowest curve at X = 0, are shown 
for reference. 



increasing the N graphene widths always results (at the same 
value of JIn) in greater superconductivity depletion in the S 
regions near the interfaces due to leakage into the N layer 
The smaller the Fermi level mismatch, the greater this effect. 
For more confining N regions (smaller Z?7v), the pair corre- 
lations decay in in the normal layer over a smaller width and 
thus the two superconductor portions of the sample are more 
strongly coupled. 

To better illustrate the slow decay of of the amplitude F{X) 
in the normal graphene region we display in Fig. 2, results for 
this normalized quantity obtained for a much larger system 
with Dn, Ds > Co- We show there also results for a broader 
range of doping levels in N. The top set of panels shows a 
global view of the correlations in the S|N|S structure, similar 
to that shown in Fig. 1, while in the bottom set of panels are 
closeups of the normal graphene region near the interface. The 
electron doped (left panels) cases are nearly identical to the 
hole doped (right panels) ones, except at smaller mismatch. 
The bottom panels allow for a more detailed examination of 
the behavior near the interface. We note that, if the mismatch 
is not large, penetration of the Cooper pairs over a distance 
clearly much larger than the correlation length occurs, and that 



depletion in the S regions occurs also in the same scale. On 
the other hand, one can see in both this and the previous fig- 
ure that the transition between the depleted S region and the 
weakly proximity-influenced N region is very abrupt: there 
are in effect two length scales, one related to the depletion 
and penetration, which can be rather longer than ^o, and an- 
other, very short scale, over which the small values of F{X) 
in N transition sharply to the depleted, but much larger, values 
in S. This is in contrast to what occurs in the standard proxim- 
ity effect in ordinary bilayer materials, which is characterized 
by a single length scale. 



Before discussing the local density of states for S |N|S struc- 
tures, it is illuminating to first investigate the DOS and char- 
acteristic energies for pure graphene nanolayers, which we do 
by setting Dg = 0. In the absence of other materials, and 
hence also of proximity effects, the DOS in this case is es- 
sentially independent of position, and thus it is appropriate to 
spatially average Eq. (3) over the entire sample, which (after 
using the normalization condition for the quasiparticle ampli- 
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FIG. 3: (Color online). Local density of states (DOS) for a normal, finite width, graphene layer that has (a) zero doping (/Iat = 0), and 
(c) moderate doping (/liv = 0.2). The normalizations of DOS and energies are chosen (see text) so that for an infinite width layer the plots 
would be straight lines of slope ±1. Four different widths are illustrated, in order of decreasing height of the main peaks: Dn ~ 20 (black). 
Dm = 100 (blue). Dm ~ 300 (red), and Dm = 600 (cyan). Each bottom panel is a magnification of the results above it, over a naiTower 
energy range. The main peaks are related to the Thouless scale, see Eq. (22) and discussion below it. 



tudes) yields a simplified normalized DOS, 



No ins) 



2Td ^ 



dk. 



cosh 



2T 



cosh 



2T 



(21) 



where T is the temperature, and iVo(e) is introduced below 
Eq. (13). In this case, with only normal material present, 
it must be understood that A^o(/^s) is just an arbitrary, but 
convenient, normalization. Note that after setting Ds = 
and hence A = 0, no iteration for self consistency is needed 
and only the eigenvalue spectra needs to be determined when 
performing the diagonalization of the matrix in Eq. (4). To 
achieve the required energy resolution for the results to follow, 
integrals over ky are numerically evaluated by transforming 
them into a sum over 5000 transverse modes. Also, in order 
to better discern the relevant DOS features, we consider the 
low temperature limit (see below). For finite width graphene 
sheets, the coherent superposition of standing waves deter- 
mines the Thouless energy scale, 



dN 



(22) 



The Thouless energy scale reveals itself in the form of peaks 
in the quasiparticle spectra that, in this simple geometry re- 
peat at odd integer multiples of Ec — ttEt- We have then 
Ec/ ^J'S = tt/Djv- The /is and kps act here as convenient ar- 
bitrary normalizations. These peaks are superimposed on the 
straight lines that would represent the DOS in the Djsi oo 
limit. If one normalizes, as we do, the energies in terms 
of and the DOS as explained above, then the slope of 
these straight lines would be ±1. Our results are shown in 
Fig. 3. The top panels ((a) and (c)) show the DOS (calcu- 
lated from and normalized as in Eq. (21)) over a broad en- 
ergy range. Panel (a) is for undoped graphene {JIn =0) and 
panel (c) corresponds to a relative doping of Jin = 0.2. For 
these two cases, there are four curves shown that correspond 
to four different graphene widths (see caption). The Thou- 
less peaks are at their predicted positions. Their magnitude 
tends to decrease as the width increases, with the results for 
largest width, Dn ~ 600, approaching the signature linear 
dispersion for bulk graphene. With the introduction of doping 
(panel (c)) the results are shifted, in normalized energy units, 
by —JiN from the Dirac point, resulting in the shifting of the 
gap away from zero energy for all widths shown. To more 
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FIG. 4: (Color online). Local density of states (DOS) (normalized as explained in the text) for an SjN|S system with highly doped S layers 
(JIn = 0). Results are shown for the S region (panels (a) and (b)) and in the N region ((c) and (d)). In (a) and (c), we have Dn = Ho/5 
(solid (black) curves) and Djv = Ho ((blue) dashed curves), while in (b) and (d) the N layers are larger: Dn = 3Ho (solid (red) curves) and 
Dn ~ 6Ho (dashed (green) curves). 



clearly discern the Thouless peaks, the bottom panels, which 
correspond to the same parameter values as the top ones, illus- 
trate the DOS over a smaller energy range. In Figs. 3 (b) and 
(d), the Thouless peaks in the DOS are clearly seen to occur 
at energies that coincide with the expression discussed below 
Eq. (22): for Dn — 300, the first peaks arise at energies corre- 
sponding to \Ec/fis\ = 7r/300 w 0.01, while for Dn = 600, 
we have, \Ec/fis\ = 7i"/600 « 0.005. In Fig. 3(b), the curves 
representing the smaller widths, and correspondingly larger 
Et, are absent since (as panel (a) shows) they emerge beyond 
the given energy window. These results illustrate also the pre- 
cision and reUablity of our methods. 



We now return to the S | N | S trilayer and investigate the roles 
that both the Thouless and superconducting energy scales play 
by considering the local DOS of a S |N|S nanostructure in both 
the S and N regions. After inserting the quasiparticle expan- 
sions found in Eq. (5), the general expression for the DOS in 



Eq. (3) can be rewritten as. 



N{x,e) Us 



dky [(1^ Sn.q sin{kqx) 



+ \''^^tn,qSill{kqX)\ j COsh 



2T 



+ ( X! sin(fcga;) + ^ 2„^g sin(fc,a;) 



n.q 

X cosh^^ 



n.q 



2T 



(23) 



In calculating the DOS for the S|N|S cases, we take the eigen- 
vectors and eigenenergies, self-consistently calculated as ex- 
plained above, and insert them into Eq. (23). When ujc < 1, 
the case considered here, the relationship between the bulk 
transition temperature Tq and Aq is found using Eq. (14), to- 
gether with,^^ 



tanh[e/(2ro)], 



(24) 



which results in the weak coupling limit in the BCS relation^^ 
Ao = (7r/7£;)To, with je being the Euler constant. We take 
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FIG. 5: (Color online). Local density of states (plotted and normalized as in Fig. 4) for an S|N|S system with moderately doped S layers 
{JiN = 0.2). The panel arrangement, curve (color and) structure, and all other parameters are the same as that in Fig. 4. See text for discussion 
and comparison 



the low temperature limit T /Tq k, 0.016. This is the same 
temperature as in the previous plots. The energy-resolved 
DOS is then determined at two locations: one at the middle 
of one of the superconducting regions, and the other at the 
center of the sample (normal region). In these plots, we will 
normalize the energy (measured as usual from the chemical 
potential) by Aq, and the plotted DOS (as before) by Nq{^s)- 
Thus, if our plots were performed for an infinite normal sam- 
ple, they would of course still be straight lines but the slope 
would now be ±Aq/iis, which is the same as l/Sg via the 
relations mentioned in Sec. 11. 

There are now two energy scales to consider One is, as 
ordinarily in all superconductors, the bulk gap Aq. In ad- 
dition, because our system is two dimensional, has a linear, 
massless dispersion relation, and is finite in the x direction, 
the DOS is also, as we have seen, drastically affected by the 
Thouless energy.-*^ The quasiparticles in the previously dis- 
cussed graphene nanostrip were confined solely by the two 
outer boundaries. Now due to the intrinsically superconduct- 
ing regions, there is also possible partial confinement by the 
self consistent pair potential A(a;), which due to Andreev 
scattering events (normal or specular), leads to a modifica- 
tion of the relevant energy scales associated with the peaks in 
the quasiparticle spectra. Thus, the two scales interact. For 



large N graphene widths (cIn ^ ^q), the normalized energy 
spectrum has the gap set primarily by the Thouless character- 
istic energy, and a peak structure'** at multiples of Et- In the 
S|N|S geometry the relevant scale is now"*^ Ec = {tt/2)Et- 
It follows then that £;c/Ao = {ttEq)/ {2Dn)- The interplay 
between this ratio and that of Aq to fis, can for a given en- 
ergy range, result in many additional resonance peaks. In the 
non self-consistent treatment, these peaks occur exactly at in- 
teger multiples of i?T,"* but as will be seen below, self con- 
sistently accounting for proximity effects in our finite sized 
system can modify this picture. Hence, as a consequence of 
the existence of the Thouless scale, the behavior of the local 
DOS in a S|N|S heterostructure is strongly dependent on the 
size of each region. If however the normal graphene channel 
is much narrower than ^o, the lowest energy scale is Aq. 

To illustrate these issues, results for Ds = 150 = I.SSq 
and various values of are exhibited in Fig. 4. These in- 
clude cases where (In < £,o and cases where it is larger, thus 
demonstrating the relative relevance of both the Aq and Et 
energy scales in different situations. Thus, consider first the 
bottom panels (c) and (d), in this Fig. 4, where results for the 
N region are plotted. In panel (d) we have = 3^o ((red) 
solid curve) and = 6^o ((green) dashed curve) so that the 
influence of the S portions of the sample is, while as we have 
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seen not negligible, weak. One clearly sees that the results can 
still be described by a straight line of slope l/Sg = 0.01 on 
which there are superimposed peaks and oscillations related 
to both the Ao and the Et scales. For the parameter values in 
this panel we have (see discussion above) that when e/Ec = 1 
then the normalized energy e/Ao is about 1/2 for (In = 3^o 
and half that for the other case shown. One sees indeed this 
behavior in this panel (d). On the other hand, in panel (c), 
where results for smaller values of Dtv/Sq are shown, the in- 
fluence of the Thouless scale is very weak for = Co ((blue) 
dashed curve) when the two energy scales roughly coincide) 
and nearly nonexistent when < £,q ((black) solid curve). 

The results in Fig. 4 for the local DOS in the S region are 
very different. These are shown in panels (a) and (b) for the 
same values of and the same plotting conventions as for 
panels (c) and (d) respectively. We now see a very clear energy 
gap close to the bulk value. The influence of the Thouless en- 
ergy is reduced now to some weak additional peaks at higher 
energy. We can see that at the energy scales shown the effect 
of the bulk normal state linear DOS is not visible, although of 
course this is an artifact arising from the energy range plotted 
and the increasing behavior reappears eventually at larger val- 
ues of e S> Aq. Even though /Iat = 0, particle-hole symmetry 
breaks down in the N region. 

The local DOS is very dependent on the doping level. 
To show this we display in Fig. 5 results for the DOS at 
piAT = 0.2. In this case, the DOS for bulk {Dn ^ oo) nor- 
mal graphene (zero Thouless energy and zero pair potential) 
is still (see Fig. 3) a straight line but with the origin shifted. 
For our parameter values, and indeed for any reasonable pa- 
rameter values in our context, this origin is shifted out of the 
horizontal scale in the energy ranges of order Aq shown in 
this figure. The four panels in the figure are arranged exactly 
as those in Fig. 4 and correspond (both the panel arrangement 
and the (color or) structure of the curves) to exactly the same 
cases. One can see that when the doping amount changes, this 
DOS becomes more complicated because of changes in the 
quasiparticle bound states shifting with the Fermi level. Thus, 
the results for large Dn in the N region (panel (d)) show now 
only a faint trace of any gap, either superconducting or Thou- 
less: the DOS is nearly linear at small energies, possessing a 
V shape at the Dirac point. This subgap structure is a mod- 
ification to the traditional Andreev bound states'"* that arise 
in the spectrum of conventional superconductor-normal metal 
systems. In panel (c), when the djv is comparable to the cor- 
relation length, the "V" behavior still persists (dashed (blue) 
curve) but it is completely gone, and replaced by a gap, when 
the thickness is below (solid (black) curve). The panels (a) 
and (b), corresponding to the S region, are less strikingly dif- 
ferent from the corresponding ones in Fig. 4 but they do show 
an intriguing additional structure in the gap region. 

In view of the strong effect, as evidenced in the compari- 
son of Figs. 4 and 5, of JIn on the gap structure in the DOS 
it is interesting to further examine in a more direct way the 
induced gap in the quasiparticle spectrum. This we do by ex- 
tracting from our numerical results the eigenvalue from the 
self consistent spectra obtained from Eq. (6) for which e„ as 
measured from the chemical potential is minimum. We call 
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FIG. 6: (Color online). The excitation gap, -Egap as defined in the 
text, as a function of the relative doping parameter, JIm- Four dif- 
ferent normal graphene widths are considered: Dn = 20, 50, 100, 
and 300 ((red) circles, (blue) squares, (green) diamonds, and (cyan) 
triangles respectively. Lines are straight segments joining points. 



this quantity the excitation gap, and denote it by Sgap, which 
is generally determined by longitudinally directed (along x) 
trajectories corresponding to small ky. The results are shown 
in Fig. 6 where we show the evolution of Sgap normalized by 
Aq (so that the quantity plotted is non-negative and less than 
unity) as a function of 'jlpf. Results for four different values 
of Dn are shown, encompassing values both above and be- 
low Sq: i^AT = 20 = 0.2So (circles), Dn ^ 50, (squares), 
Dn = 100 (diamonds), and Dn = 300 (triangles). In all 
cases Ds = 150. The range of JIn and values of Dn tiiat 
result in a gap are of course consistent with the DOS results 
above. The results shown illustrate that structures including 
narrower normal graphene layers possess energy gaps that are 
much more robust to changes in N layer doping. The contrac- 
tion of the gap with increasing Dn is qualitatively similar to 
what is observed in conventional three dimensional systems, 
but, as mentioned above, the structure of the gap amplitude 
and of the DOS is very different. 

Up to this point, we have considered the low temperature 
limit. It is of interest both experimentally and theoretically to 
now turn our attention to the calculation of the critical tem- 
perature Tc of the S|N|S structures and its dependence on 
doping levels and geometrical parameters. This quantity is 
calculated using the efficient eigenvalue method described by 
Eq. (16) and the discussion below it. Results are presented in 
terms of the ratio Tc/Tq and displayed in Fig. 7. In the left 
panel, results are given as a function of relative doping level 
pi AT for two values of Ds'. Ds = Sq = 100 ((blue) squares) 
and Ds = So/2 ((red) circles). We keep Dn — con- 
stant in this figure. We see than increasing IJInI to moderate 
values, that is, decreasing the Fermi level mismatch, leads to 
(Fig. 7(a)) a reduction in via the corresponding increase in 
the interlayer coupling. This effect is more pronounced for 
thin S layers, where Tc can vary with JIn in a nontrivial fash- 
ion. It is remarkable, however, that Tc remains rather high 
even when Ds is smaller than the correlation length. This is 
is stark contrast to ordinary three dimensional S |N|S systems, 
where Tc drops much more rapidly for small thicknesses satis- 
fying ds ^ £,0- There is also a clear asymmetry in the critical 
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FIG. 7: (Color online). The critical temperature Tc (normalized by 
To) for a S|N|S system as a function of (a) ptjv, and (b) the S width, 
Ds- In (a) the doping dependence is shown for two values of Ds ~ 
50 ((red) circles) and 100 ((blue) squares). In (b), results are shown 
as a function of Ds for three doping levels: JIn = ((red) circles), 
0.2 ((blue) squares) and 0.4 ((green) diamonds). The doping level 
in these cases has a dramatic effect on Tc for small Ds, but is less 
detrimental for larger Ds, where the curves tend to coalesce in the 
limit of bulk S widths. In both (a) and (b), the normal graphene layer 
has Dn = 100. 



temperature as a function of doping, where for a given mag- 
nitude \]j.n\, electron doping more strongly reduces Tc- In the 
right panel, (Fig. 7(b)) we illustrate that varying the width of 
the superconductors has a considerably greater impact on Tc 
for moderate values of JIn than when the mismatch is large. In 
the latter case the Ds dependence remains weak as long as Ds 
is still comparable to Sq, however (and consistent with panel 

(a) ) the superconducting regions that have widths a fraction of 
the coherence length reveal the richest behavior. Continuing 
to reduce Ds beyond some critical value, of course results in 
the graphene system eventually becoming nonsuperconduct- 
ing, as Cooper pair formation is inhibited. 

Results such as those shown in Fig. 7(a) imply that at a 
fixed temperature, variations in the doping parameter, Jlisi, can 
lead to a S|N|S system transitioning from a superconducting 
state to normal one and vice versa. Since graphene doping 
can be effecttively tuned via application of an external electric 
field,'"^ this may offer possibilities as a carbon-based S|N|S 
switch for supercurrent flow. This question is, therefore, wor- 
thy of further discussion. We thus expand on this point by 
showing in Fig. 8 the normalized pair ampUtude F{X) plot- 
ted as a function of X for several positive values of ptjv. Each 
panel is at a different fixed temperature and the geometrical 
parameters, Ds = So/2 and Ds ~ So are chosen to corre- 
late respectively with the (red) circled and (blue) squared data 
of Fig. 7(a). The two representative temperatures that we in- 
vestigate are T = O.STTq (panel (a)), and T = 0.92To (panel 

(b) ). The graphene region that is intrinsically nonsupercon- 
ducting has a width in both cases corresponding to Dn = Eq. 
One can see in Fig. 7(a) that for the smaller Ds ~ So/2, the 
temperature T = O.STTo corresponds to Tc near JIn ^ 0.35, 
and for Ds = Sq, the temperature T = 0.92ro results in Tc 
near JIn ~ 0.4. The regions for which positive JIn is smaller 
being superconducting (the corresponding negative JIn differ 




FIG. 8: (Color online). The normalized pair amplitude as a function 
of position. In (a), the temperature is set at T = O.STTo, and the 
normal and superconductor widths are Djv = 100 and Ds = 50 
respectively. In (b), we have T = 0.92ro, while Dn = 100 and 
Ds ~ 100. In both cases, the doping parameter, ptjv, is varied from 
to 0.4 in increments of 0.1. The aiTow depicts the direction of 
increasing JIm which eventually leads to the vanishing of the pairing 
correlations. 



slightly due to the electron-hole doping asymmetry). This is 
more clearly seen in Fig. 8, where as Jin is increased, the pair 
amplitude is seen to decrease before plummeting abruptly to 
zero as JIn reaches its critical value: near JIn —?■ 0.35 in panel 
(a) or near JIn — >■ 0.4 for panel (b). Thus, if this transition can 
be manipulated via electric fields, abrupt switching will result. 



4. CONCLUSIONS 

We have studied in this paper the proximity effects that 
occur in clean, doped and undoped, graphene-based S|N|S 
trilayers. We have created and implemented a fully self- 
consistent procedure to calculate the electron and hole wave- 
functions and energy spectrum of the system, from which we 
have extracted the pair amplitude and the local DOS. We also 
developed a semianalytical and computationally efficient lin- 
earized method that can calculate the transition temperature, 
Tc, of the system. 
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We have found that the behavior of the pair amplitude near 
the interfaces (the directly observed proximity effect) depends 
strongly on the relative doping levels of the S and the N por- 
tions, and that the pair amplitude is described by two differ- 
ent length scales. One length scale is related to penetration 
of the superconducting correlations, and is long ranged (rel- 
ative to ^o) and the other scale is short ranged, and corre- 
lates to Cooper pair leakage from the S regions near the in- 
terfaces. We illustrated that if the normal graphene layer is 
weakly doped, specular Andreev reflection can lead to super- 
conducting correlations penetrating into the normal graphene 
region. The local DOS exhibits a number of striking features, 
arising from the interplay between the superconducting and 
the Thouless energy scales. This interplay depends of course 
on geometry, where the two energy scales overlap when the 
graphene layer widths are on the same order as ^o- For our 
larger structures (with widths exceeding ^o)> undoped normal 
regions revealed resonant peaks and energy gaps at character- 
istic energies proportional to the Thouless energy scale Et- 
By moderately doping the N region, there was an emergence 
of Andreev bound states in the S regions and a destruction 
of the energy gap. The smaller S|N|S structures (with widths 
smaller than ^o) revealed energy gaps that are linked mainly 
to the Aq scale, and are more robust to doping. 

We also developed a general microscopic method for calcu- 
lating Tc for S|N|S nanostructures, by linearizing the DBdG 
equations and the self consistency condition. We found that 
for small S layer widths, decreasing the Fermi level mismatch 
leads to a nontrivial reduction in Tc- The critical temperature 



also exhibited a clear asymmetry as a function of doping, and 
typically electron doping had a greater impact on reducing Tc- 
Thus if doping is to be modified by an electric field, ''^^ the 
polarity"*" of the field can have an important effect on the crit- 
ical temperature. The study of Tc revealed reentrant behavior 
as a function of doping. These behaviors may lead to switch- 
ing phenomena as a function of applied electric field, and thus 
depending on the bias, superconductivity can be turned on or 
off. The effectiveness of graphene as a low temperature field 
effect device therefore depends in large part by the proxim- 
ity effects, which can only be accounted for within a self- 
consistent framework. This work represents the first step, a 
proof of principle, as to the use of our self consistent methods 
in graphene. Other issues, such as those related to ferromag- 
netically doped graphene in contact with a superconductor re- 
gion, can also be examined using the same techniques. We 
expect that many aspects of the ever intriguing behavior of 
graphene-based heterostructures will be illuminated via appli- 
cation of these methods. 
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